library(data.table)
library(mefa)
library(doParallel)
library(emayili)
library(magrittr)
library(ggplot2)
library(nleqslv)



source('helper_functions.r')

runBootstraps <- function(estimation) {
  
  N = 100
  for(NN in 1:N) {
    runCounterfactualBootstrap(estimation,paste0('bootstrap_',as.numeric(Sys.time())))
  }
  
}


combineBootstraps <- function() {
  
  allFiles <- dir('counterfactual_bootstrap/all_estimation_24/',full.names = T)
  
  for(file.index in 1:length(allFiles)) {
    load(allFiles[file.index])
    cr[,run := file.index]
    qe[,run := file.index]
    lim[,run := file.index]
    if(file.index == 1) {
      all.runs.cr <- cr
      all.runs.qe <- qe
      all.runs.lim <- lim
    } else {
      all.runs.cr <- rbind(all.runs.cr,cr)
      all.runs.qe <- rbind(all.runs.qe,qe)
      all.runs.lim <- rbind(all.runs.lim,lim)
    }
  }
  
  # Labels
  labels <- data.table(variable = c('3pct','45pct','baseline','75pct','9pct','12pct'),req = c(.03,.045,.06,.075,.09,.12))
  
  # Total volume
  ff <- melt(all.runs.cr[variable == 'total.volume',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_line(aes(x=req,y=med)) + geom_ribbon(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    coord_cartesian(ylim = c(1250,2000)) + 
    xlab('Capital requirement') + 
    theme_bw() + 
    ylab('Lending volume') + 
    scale_x_continuous(breaks = labels$req)
  ggsave('counterfactual_bootstrap/total_lending.png',height=4,width=6,units = 'in')
  
  # Balance sheet financing share
  ff <- melt(all.runs.cr[variable == 'balance.sheet.financing.share',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_line(aes(x=req,y=med)) + geom_ribbon(aes(x=req,ymin = p025,ymax=p975),alpha = .25)  + 
    coord_cartesian(ylim = c(0,1)) + 
    xlab('Capital requirement') + 
    theme_bw() + 
    ylab('Balance sheet financing share') + 
    scale_x_continuous(breaks = labels$req)
  ggsave('counterfactual_bootstrap/balance_sheet_financing_share.png',height=4,width=6,units = 'in')
  
  # Shadow bank market share
  ff <- melt(all.runs.cr[variable == 'shadow.bank.market.share',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_line(aes(x=req,y=med)) + geom_ribbon(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    coord_cartesian(ylim = c(0,1)) + 
    xlab('Capital requirement') + 
    theme_bw() + 
    ylab('Shadow bank share') + 
    scale_x_continuous(breaks = labels$req)
  ggsave('counterfactual_bootstrap/shadow_bank_share.png',height=4,width=6,units = 'in')
  
  # Rate spread
  ff <- melt(all.runs.cr[variable == 'rate.spread',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_line(aes(x=req,y=med)) + geom_ribbon(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    #coord_cartesian(ylim = c(0,1)) + 
    xlab('Capital requirement') + 
    theme_bw() + 
    ylab('Jumbo - conforming rate spread') + 
    scale_x_continuous(breaks = labels$req)
  ggsave('counterfactual_bootstrap/rate_spread.png',height=4,width=6,units = 'in')
  
  
  
  # Labels
  labels <- data.table(variable = c('m100','m25','m10','baseline','p10','p25','p100'),req = c(-100,-25,-10,0,10,25,100))
  
  # Total volume
  ff <- melt(all.runs.qe[variable == 'total.volume',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_line(aes(x=req,y=med)) + geom_ribbon(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    xlab('GSE financing costs (bps)') + 
    theme_bw() + 
    ylab('Lending volume') + 
    scale_x_continuous(breaks = labels$req)
  ggsave('counterfactual_bootstrap/total_lending_qe.png',height=4,width=6,units = 'in')
  
  # Balance sheet financing share
  ff <- melt(all.runs.qe[variable == 'balance.sheet.financing.share',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_line(aes(x=req,y=med)) + geom_ribbon(aes(x=req,ymin = p025,ymax=p975),alpha = .25)  + 
    xlab('GSE financing costs (bps)') + 
    theme_bw() + 
    ylab('Balance sheet financing share') + 
    scale_x_continuous(breaks = labels$req)
  ggsave('counterfactual_bootstrap/balance_sheet_financing_share_qe.png',height=4,width=6,units = 'in')
  
  # Shadow bank market share
  ff <- melt(all.runs.qe[variable == 'shadow.bank.market.share',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_line(aes(x=req,y=med)) + geom_ribbon(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    xlab('GSE financing costs (bps)') + 
    theme_bw() + 
    ylab('Shadow bank share') + 
    scale_x_continuous(breaks = labels$req)
  ggsave('counterfactual_bootstrap/shadow_bank_share_qe.png',height=4,width=6,units = 'in')
  
  # Rate spread
  ff <- melt(all.runs.qe[variable == 'rate.spread',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_line(aes(x=req,y=med)) + geom_ribbon(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    xlab('GSE financing costs (bps)') + 
    theme_bw() + 
    ylab('Jumbo - conforming rate spread') + 
    scale_x_continuous(breaks = labels$req)
  ggsave('counterfactual_bootstrap/rate_spread_qe.png',height=4,width=6,units = 'in')
  
  
  # Labels
  labels <- data.table(variable = c('m25','baseline','p25','none','all_low','all_high'),req = c('1. -25%','2. Baseline','3. +25%','4. None','5. All 417k','6. All 625k'))
  
  # Total volume
  ff <- melt(all.runs.lim[variable == 'total.volume',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_point(aes(x=req,y=med),stat = 'identity') + geom_errorbar(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    xlab('Conforming loan limit') + 
    theme_bw() + 
    ylab('Lending volume') + 
    scale_x_discrete(breaks = labels$req)
  ggsave('counterfactual_bootstrap/total_lending_lim.png',height=4,width=6,units = 'in')
  
  # Balance sheet financing share
  ff <- melt(all.runs.lim[variable == 'balance.sheet.financing.share',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_point(aes(x=req,y=med),stat = 'identity') + geom_errorbar(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    xlab('Conforming loan limit') + 
    theme_bw() + 
    ylab('Balance sheet financing share') + 
    scale_x_discrete(breaks = labels$req)
  ggsave('counterfactual_bootstrap/balance_sheet_financing_share_lim.png',height=4,width=6,units = 'in')
  
  # Shadow bank market share
  ff <- melt(all.runs.lim[variable == 'shadow.bank.market.share',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_point(aes(x=req,y=med),stat = 'identity') + geom_errorbar(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    xlab('Conforming loan limit') + 
    theme_bw() + 
    ylab('Shadow bank share') + 
    scale_x_discrete(breaks = labels$req)
  ggsave('counterfactual_bootstrap/shadow_bank_share_lim.png',height=4,width=6,units = 'in')
  
  # Rate spread
  ff <- melt(all.runs.lim[variable == 'rate.spread',!c('variable')],id.vars = 'run')
  rr <- ff[,j=list(p025 = quantile(value,.025,na.rm=T),med = median(value,na.rm=T),p975 = quantile(value,.975,na.rm=T)),by='variable']
  rr <- merge(rr,labels,by='variable')
  ggplot(rr) + geom_point(aes(x=req,y=med),stat = 'identity') + geom_errorbar(aes(x=req,ymin = p025,ymax=p975),alpha = .25) + 
    xlab('Conforming loan limit') + 
    theme_bw() + 
    ylab('Jumbo - conforming rate spread') + 
    scale_x_discrete(breaks = labels$req)
  ggsave('counterfactual_bootstrap/rate_spread_lim.png',height=4,width=6,units = 'in')
  
}


runCounterfactualBootstrap <- function(estimation,filename) {
  
  ## Get the bootstrapped parameters ready
  source(paste0('estimation_configs/',estimation,'.r'))
  load(paste0('estimation_results/',pars$filename,'.rsave'))
  
  # Get all the deltas and market shares
  
  # 0. Get estimation pars
  pars <- estimation_data$pars
  
  # 1.  Deal parameters
  new_mx <- param_vec_to_matrix(estimation_data$solution$x)
  A <- new_mx$A
  B <- new_mx$B
  S <- new_mx$S
  
  linear <- estimation_data$solution$linear.model
  
  # Standoard errors on parameters
  se_alpha = 0.27
  se_beta = 0.14
  se_gamma = 0.13
  se_log_f = 0.03
  
  se_alpha_income = 0.0017
  se_alpha_price  = 0.0014
  se_beta_income  = 0.0019
  se_beta_price   = 0.0018
  se_gamma_income = 0.0004
  se_gamma_price  = 0.0005
  se_log_f_income = 0.0037
  se_log_f_price  = 0.0035
  
  se_alpha_sig = 0.0005
  se_beta_sig  = 0.0001
  se_gamma_sig = 0.0001
  se_log_f_sig = 0.0112
  
  A['beta','A']  <- pars$A['beta','A']  + rnorm(1,0,se_beta)
  A['log_f','A'] <- pars$A['log_f','A'] + rnorm(1,0,se_log_f)
  
  B['alpha','log_income']      <- pars$B['alpha','log_income']      + rnorm(1,0,se_alpha_income)
  B['alpha','log_house_price'] <- pars$B['alpha','log_house_price'] + rnorm(1,0,se_alpha_price)
  B['beta','log_income']       <- pars$B['beta','log_income']       + rnorm(1,0,se_beta_income)
  B['beta','log_house_price']  <- pars$B['beta','log_house_price']  + rnorm(1,0,se_beta_price)
  B['gamma','log_income']      <- pars$B['gamma','log_income']      + rnorm(1,0,se_gamma_income)
  B['gamma','log_house_price'] <- pars$B['gamma','log_house_price'] + rnorm(1,0,se_gamma_price)
  B['log_f','log_income']      <- pars$B['log_f','log_income']      + rnorm(1,0,se_log_f_income)
  B['log_f','log_house_price'] <- pars$B['log_f','log_house_price'] + rnorm(1,0,se_log_f_price)
  
  S['alpha','alpha'] <- pars$S['alpha','alpha'] + rnorm(1,0,se_alpha_sig)
  S['beta','beta']   <- pars$S['beta','beta']   + rnorm(1,0,se_beta_sig)
  S['gamma','gamma'] <- pars$S['gamma','gamma'] + rnorm(1,0,se_gamma_sig)
  S['log_f','log_f'] <- pars$S['log_f','log_f'] + rnorm(1,0,se_log_f_sig)
  
  linear$coefficients = linear$coefficients + coef(summary(linear))[, "Std. Error"] * rnorm(length(linear$coefficients),0,1)
  
  ## Supply estimation bootstrapping
  load(paste0('supply_estimation_results/',estimation,'.rsave'))
  
  supply_linear = result$mc$linear
  supply_nonlinear = table.D
  supply_year <- table.A[,c('year','value')]
  
  supply_linear$coefficients = supply_linear$coefficients + coef(summary(supply_linear))[, "Std. Error"] * rnorm(length(supply_linear$coefficients),0,1)
  
  supply_nonlinear$value[[1]] <- supply_nonlinear$value[[1]] + supply_nonlinear$se[[1]] * rnorm(1,0,1)
  supply_nonlinear$value[[2]] <- supply_nonlinear$value[[2]] + supply_nonlinear$se[[2]] * rnorm(1,0,1)
  supply_nonlinear$value[[3]] <- supply_nonlinear$value[[3]] + supply_nonlinear$se[[3]] * rnorm(1,0,1)
  
  supply_year$value <- supply_year$value + table.A$se * rnorm(nrow(supply_year),0,1)
  names(supply_year) <- c('YEAR','MC_LIN_YEAR')
  
  
  ENTRY = T
  BANK_EXIT = F
  
  # Get baseline scale
  baseline <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  SCALE = 1710 / baseline[variable == 'total.volume']$baseline
  
  # Conforming loan limits
  cr.30 <- runCounterfactual(estimation = estimation,CR = 0.030,LIM = NULL,RGSE = NULL,NAME = '3pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  cr.45 <- runCounterfactual(estimation = estimation,CR = 0.045,LIM = NULL,RGSE = NULL,NAME = '45pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  cr.60 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  cr.75 <- runCounterfactual(estimation = estimation,CR = 0.075,LIM = NULL,RGSE = NULL,NAME = '75pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  cr.90 <- runCounterfactual(estimation = estimation,CR = 0.090,LIM = NULL,RGSE = NULL,NAME = '9pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  cr.12 <- runCounterfactual(estimation = estimation,CR = 0.120,LIM = NULL,RGSE = NULL,NAME = '12pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  cr = merge(cr.30,merge(cr.45,merge(cr.60,merge(cr.75,merge(cr.90,cr.12,by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  # define rel columns
  rel.columns = c('rate.conforming','rate.jumbo','rate.spread','lender.profit','bank.profit','sb.profit','overall.surplus','individual.surplus','top.market.surplus','bottom.market.surplus','individual.surplus.top','individual.surplus.low')
  cr[variable %in% rel.columns,`3pct` := `3pct` - baseline ]
  cr[variable %in% rel.columns,`45pct` := `45pct` - baseline ]
  cr[variable %in% rel.columns,`75pct` := `75pct` - baseline ]
  cr[variable %in% rel.columns,`9pct` := `9pct` - baseline ]
  cr[variable %in% rel.columns,`12pct` := `12pct` - baseline ]
  cr[variable %in% rel.columns, baseline := 0 ]
  
  # QE
  qe.m.100 <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = -1.00,NAME = 'm100',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  qe.m.25  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = -0.25,NAME = 'm25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  qe.m.10  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = -0.10,NAME = 'm10',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  qe.base  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  qe.p.10  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE =  0.10,NAME = 'p10',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  qe.p.25  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE =  0.25,NAME = 'p25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  qe.p.100 <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE =  1.00,NAME = 'p100',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  qe = merge(qe.m.100,merge(qe.m.25,merge(qe.m.10,merge(qe.base,merge(qe.p.10,merge(qe.p.25,qe.p.100,by='variable'),by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  qe[variable %in% rel.columns, m100 := m100 - baseline ]
  qe[variable %in% rel.columns, m25 := m25 - baseline ]
  qe[variable %in% rel.columns, m10 := m10 - baseline ]
  qe[variable %in% rel.columns, p10 := p10 - baseline]
  qe[variable %in% rel.columns, p25 := p25 - baseline ]
  qe[variable %in% rel.columns, p100 := p100 - baseline ]
  qe[variable %in% rel.columns, baseline := 0 ]
  
  
  # Limits
  lim.m.25   <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 0.75,RGSE = NULL,NAME = 'm25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  lim.base   <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  lim.p.25   <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 1.25,RGSE = NULL,NAME = 'p25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  lim.p.none <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 5   ,RGSE = NULL,NAME = 'none',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  lim.lower  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 'all_low',RGSE = NULL,NAME = 'all_low',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  lim.upper  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 'all_high',RGSE = NULL,NAME = 'all_high',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,override_A = A,override_B = B,override_S = S,override_linear_demand = linear,override_linear_year = supply_year,override_linear_supply = supply_linear,override_nonlinear_supply = supply_nonlinear)
  lim = merge(lim.m.25,merge(lim.base,merge(lim.p.25,merge(lim.p.none,merge(lim.lower,lim.upper,by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  lim[variable %in% rel.columns, m25 := m25 - baseline ]
  lim[variable %in% rel.columns, p25 := p25 - baseline ]
  lim[variable %in% rel.columns, all_low := all_low - baseline ]
  lim[variable %in% rel.columns, none := none - baseline ]
  lim[variable %in% rel.columns, all_high := all_high - baseline]
  lim[variable %in% rel.columns, baseline := 0 ]
  
  
  save(list = c('cr','qe','lim'),file = paste0('counterfactual_bootstrap/all_',estimation,'/',filename,'.rsave'))  
}


runCounterfactual <- function(estimation = 'estimation_16',CR = NULL,LIM = NULL,RGSE = NULL,EQUITY_ISSUANCE_RATE = 999,JUMBO_SEC_RATE = 999, NAME = 'baseline',scale = NULL,entry = F,ALLOW_BANK_EXIT = F,risk = 'ALL',
                              override_A = NULL,override_B = NULL,override_S = NULL,override_linear_demand = NULL,
                              override_linear_year = NULL,override_linear_supply = NULL,override_nonlinear_supply = NULL) {
  
  
  cf.data <- prepCounterfactualData(estimation,CR,LIM,RGSE,EQUITY_ISSUANCE_RATE = EQUITY_ISSUANCE_RATE,JUMBO_SEC_RATE = JUMBO_SEC_RATE,ALLOW_BANK_EXIT = ALLOW_BANK_EXIT,
                                    override_A = override_A,override_B = override_B,override_S = override_S,override_linear_demand = override_linear_demand,
                                    override_linear_year = override_linear_year,override_linear_supply = override_linear_supply,override_nonlinear_supply = override_nonlinear_supply)
  if(risk == 'LOW') {
    cf.data$market.data <- cf.data$market.data[FICO == 'l']
  } else if(risk == 'HIGH') {
    cf.data$market.data <- cf.data$market.data[FICO == 'h']
  }
  
  markets.data <- cf.data$market.data[YEAR == 2017]
  markets.data[,max.xi := max(abs(xi)),by='MARKET_ID']
  markets.data <- markets.data[max.xi < 10 ]
  markets.data[,MARKET_INCOME_BUCKET := 'middle']
  markets.data[market.log_income.mean < quantile(markets.data$market.log_income.mean,.25),MARKET_INCOME_BUCKET := 'low']
  markets.data[market.log_income.mean > quantile(markets.data$market.log_income.mean,.75),MARKET_INCOME_BUCKET := 'high']
  
  if(entry) {
    solved <- getAllEquilbriumInterestRatesWithEntry(markets.data = markets.data,linear.model = cf.data$linear.model,raw.draws = cf.data$raw.draws,demand.parameters = cf.data$demand.pars)  
  } else {
    solved <- getAllEquilbriumInterestRates(markets.data = markets.data,linear.model = cf.data$linear.model,raw.draws = cf.data$raw.draws,demand.parameters = cf.data$demand.pars)  
  }
  
  
  # Merge it with info about lenders
  ll <- merge(solved,markets.data[,c('MARKET_ID','YEAR','CBSA','PURPOSE','FICO','j','LENDER','JUMBO','HELD','MARKET_SIZE','MARKET_SHARE_DATA','MARKET_INCOME_BUCKET','ALTERNATIVE')],by=c('MARKET_ID','j'))
  
  # Get total size for adjustment
  if(!is.null(scale)) {
    ll[,SHARE_MODEL := SHARE_MODEL * scale]
  }

  weighted.mean(ll$RATE_INST - ll$MC,w=ll$MARKET_SIZE * ll$SHARE_MODEL,na.rm=T)
  
  # Stats
  volumes.m = ll[,j=list(total.volume = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE,na.rm=T)/1e9,
                       conforming.volume = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c'),na.rm=T)/1e9,
                       jumbo.volume = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'j'),na.rm=T)/1e9,
                       bank.volume  = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (LENDER == 'b'),na.rm=T)/1e9)]
  
  financing.m  = ll[,j=list(balance.sheet.financing = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * HELD,na.rm=T)/1e9,
                            balance.sheet.financing.share = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * HELD,na.rm=T) / sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE,na.rm=T),
                            conforming.balance.sheet.share = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c') * HELD,na.rm=T) / sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c'),na.rm=T),
                            shadow.bank.market.share = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (LENDER != 'b') ,na.rm=T) / sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE ,na.rm=T),
                            shadow.bank.conforming.share = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c') * (LENDER != 'b') ,na.rm=T) / sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c'),na.rm=T),
                            alternative_financing = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * ALTERNATIVE,na.rm=T)/1e9)]
  
  rates.m = ll[,j=list(rate.conforming = .01 * weighted.mean(RATE_INST,w = MARKET_SIZE * SHARE_MODEL * (JUMBO == 'c'),na.rm=T),
                       rate.jumbo = .01 * weighted.mean(RATE_INST,w = MARKET_SIZE * SHARE_MODEL * (JUMBO == 'j'),na.rm=T),
                       rate.spread = .01 * weighted.mean(RATE_INST,w = MARKET_SIZE * SHARE_MODEL * (JUMBO == 'j'),na.rm=T) -  .01 * weighted.mean(RATE_INST,w = MARKET_SIZE * SHARE_MODEL * (JUMBO == 'c'),na.rm=T))]
  
  profits.m = ll[,j=list(lender.profit = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (RATE_INST - MC)/100 * ( (1- (1+0.05)^(-10) ) / (0.05)),na.rm=T)/1e9,
                         bank.profit   = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (RATE_INST - MC)/100 * ( (1- (1+0.05)^(-10) ) / (0.05)) * (LENDER == 'b')/1e9,na.rm=T),
                         sb.profit   = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (RATE_INST - MC)/100 * ( (1- (1+0.05)^(-10) ) / (0.05)) * (LENDER != 'b')/1e9,na.rm=T),
                         mean.nfsb   = weighted.mean(N_LENDERS,w=as.numeric(LENDER == 'n'),na.rm=T))]
  
  # Surplus is in units of surplus per dollar of loan...
  # First, drop the repeats (should be equal by market)
  ll.surplus <- ll[!duplicated(MARKET_ID)]
  surplus.m = ll.surplus[,j=list(overall.surplus = sum(SURPLUS * SHARE_MODEL * MARKET_SIZE * LOAN_SIZE,na.rm=T)/1e9,
                         individual.surplus = weighted.mean(SURPLUS * SHARE_MODEL * LOAN_SIZE,w=,na.rm = T),
                         top.market.surplus = sum(SURPLUS * SHARE_MODEL * MARKET_SIZE * (MARKET_INCOME_BUCKET == 'high') * LOAN_SIZE,na.rm=T )/1e9,
                         bottom.market.surplus = sum(SURPLUS * SHARE_MODEL * MARKET_SIZE * (MARKET_INCOME_BUCKET == 'low') * LOAN_SIZE,na.rm=T )/1e9,
                         individual.surplus.top = weighted.mean(SURPLUS_HIGH * SHARE_MODEL * LOAN_SIZE,na.rm=T),
                         individual.surplus.low = weighted.mean(SURPLUS_LOW * SHARE_MODEL * LOAN_SIZE,na.rm=T))]
  
  all.columns = cbind(volumes.m,financing.m,rates.m,profits.m,surplus.m)
  result <- melt(all.columns)
  names(result) = c('variable',NAME)
  
  return(result)
}

getAllEquilbriumInterestRates <- function(markets.data,linear.model,raw.draws,demand.parameters) {
  
  ## Parallelized and non-parallelized
  all.markets <- unique(markets.data$MARKET_ID)
  
  CORES = 30
  
  ## Parallelized
  if(CORES > 1) {
    registerDoParallel(cores = CORES)
    solved.model <- foreach(ii = 1:length(all.markets), .combine = rbind) %dopar% {
        current.market.id <- all.markets[ii]
        current.market.data       <- markets.data[MARKET_ID == current.market.id]
        getEquilibriumInterestRate(current.market.data,linear.model,raw.draws,demand.parameters)
    }
  } else {  
    for(ii in 1:length(all.markets)) {
      current.market.id <- all.markets[ii]
      print(current.market.id)
      current.market.data       <- markets.data[MARKET_ID == current.market.id]
      if(ii == 1) {
        solved.model <- getEquilibriumInterestRate(current.market.data,linear.model,raw.draws,demand.parameters)
      } else {
        solved.model <- rbind(solved.model,getEquilibriumInterestRate(current.market.data,linear.model,raw.draws,demand.parameters))
      }
      
    } 
  }
  
  
  return(solved.model) 
  
}


getAllEquilbriumInterestRatesWithEntry <- function(markets.data,linear.model,raw.draws,demand.parameters) {
  
  ## Parallelized and non-parallelized
  all.markets <- unique(markets.data$MARKET_ID)
  
  CORES = 30
  
  ## Parallelized
  if(CORES > 1) {
    registerDoParallel(cores = CORES)
    solved.model <- foreach(ii = 1:length(all.markets), .combine = rbind) %dopar% {
      current.market.id <- all.markets[ii]
      current.market.data       <- markets.data[MARKET_ID == current.market.id]
      getEquilibriumInterestRateWithEntry(current.market.data,linear.model,raw.draws,demand.parameters)
    }
  } else {  
    for(ii in 1:length(all.markets)) {
      current.market.id <- all.markets[ii]
      print(current.market.id)
      current.market.data       <- markets.data[MARKET_ID == current.market.id]
      if(ii == 1) {
        solved.model <- getEquilibriumInterestRateWithEntry(current.market.data,linear.model,raw.draws,demand.parameters)
      } else {
        solved.model <- rbind(solved.model,getEquilibriumInterestRateWithEntry(current.market.data,linear.model,raw.draws,demand.parameters))
      }
      
    } 
  }
  
  
  return(solved.model) 
  
}

getEquilibriumInterestRate <- function(market.data,linear.model,raw.draws,demand.parameters) {
  
  market.data <- market.data[order(j)]
  transformed.draws <- transformRawDraws(raw.draws,market.data,demand.parameters)
  
  objective.function <- function(x,info = F) {
    
    # Assume you get 'x' in the right order
    market.data$RATE_INST <- x
    
    # Get markup
    mr <- calculateMarginalRevenuesAndCounterfactualShares(market.data = market.data,transformed.draws = transformed.draws,pars = demand.parameters,linear.model = linear.model)
    mr <- merge(mr,market.data[,c('j','RATE_INST','MC','N_LENDERS')],by='j')
    mr[,foc := MC + MARKUP_MODEL - RATE_INST] 
    val <- mr$foc
    if(info) {
      return(mr[,c('MARKET_ID','j','SHARE_MODEL','LOAN_SIZE','MARKUP_MODEL','RATE_INST','MC','foc','SURPLUS','SURPLUS_HIGH','SURPLUS_LOW','N_LENDERS'),with=F])
    } else {
      return(val)
    }
  }
  
  
  soln <- nleqslv(x = market.data$RATE_INST,fn = objective.function)
  result <- objective.function(soln$x,info = T)  
  return(result)  
}


getEquilibriumInterestRateWithEntry <- function(market.data,linear.model,raw.draws,demand.parameters) {
  
  market.data <- market.data[order(j)]
  transformed.draws <- transformRawDraws(raw.draws,market.data,demand.parameters)
  
  cost_m = 10.37
  cost_s = 7.30
  cost_N = 815
  
  # This needs tobe fixed to (1) more flexibly accomodate different parameters
  # and (2) incorporate the presense of n and f types doing jumbos...
  
  
  objective.function <- function(x,info = F) {
    
    # Assume you get 'x' in the right order
    toWork <- market.data
    
    
    toWork$RATE_INST <- x[1:nrow(market.data)]
    toWork[LENDER == 'n',N_LENDERS := x[1 + nrow(market.data)]]
    jj = toWork[LENDER == 'n']$j
    
    # Get markup
    mr <- calculateMarginalRevenuesAndCounterfactualShares(market.data = toWork,transformed.draws = transformed.draws,pars = demand.parameters,linear.model = linear.model)
    mr <- merge(mr,toWork[,c('j','RATE_INST','MC','MARKET_SIZE','N_LENDERS')],by='j')
    mr[,foc := MC + MARKUP_MODEL - RATE_INST] 
    mr[,VARIABLE_PROFIT := MARKET_SIZE * SHARE_MODEL/N_LENDERS * LOAN_SIZE * (MARKUP_MODEL / 100)/1000 ] # in k dollars, per lender]
    nfsb.var.profit = max(0,sum(mr[j %in% jj]$VARIABLE_PROFIT))
    nfsb.E_LENDERS = cost_N * pnorm(log(nfsb.var.profit),cost_m,cost_s)
    D_LENDERS =  nfsb.E_LENDERS - x[1 + nrow(market.data)]
    val <- c(mr$foc,D_LENDERS)
    if(info) {
      return(mr[,c('MARKET_ID','j','SHARE_MODEL','LOAN_SIZE','MARKUP_MODEL','RATE_INST','MC','foc','SURPLUS','SURPLUS_HIGH','SURPLUS_LOW','N_LENDERS'),with=F])
    } else {
      return(val)
    }
  }
  
  
  soln <- nleqslv(x = c(market.data$RATE_INST,market.data[LENDER == 'n']$N_LENDERS[1]),fn = objective.function)
  result <- objective.function(soln$x,info = T)  
  return(result)  
}



prepCounterfactualData <- function(estimation = 'estimation_5',CR = NULL,LIM = NULL,RGSE = NULL,JUMBO_SEC_RATE = 999, EQUITY_ISSUANCE_RATE = 999,ALLOW_BANK_EXIT = F,
                                   override_A = NULL,override_B = NULL, override_S = NULL,override_linear_demand = NULL,
                                   override_linear_year = NULL,override_linear_supply = NULL, override_nonlinear_supply = NULL) {
  
  
  
  # Load demand estimation data
  demand.estimation <- loadEstimation(estimation,override_A = override_A,override_B = override_B,override_S = override_S, override_linear_demand = override_linear_demand)
  market.data       <- demand.estimation$market.data
  linear.model      <- demand.estimation$linear.model
  demand.pars       <- demand.estimation$pars
  raw.draws         <- demand.estimation$raw.draws
  
  
  # Counterfactual conforming loan limit?
  if(!is.null(LIM)) {
    if(LIM == 'all_low') {
      market.data[, market.conforming_loan_limit := 417000]
    } else if(LIM == 'all_high') {
      market.data[, market.conforming_loan_limit := 625000]
    } else {
      market.data[, market.conforming_loan_limit := market.conforming_loan_limit * LIM]
    }
  }
  
  # Load supply estimation data
  supply.estimation <- loadSupplyEstimationForCF(estimation,override_linear_year = override_linear_year,
                                                            override_linear_supply = override_linear_supply,
                                                            override_nonlinear_supply = override_nonlinear_supply)
  estimation.data <- loadSupplyEstimationData(estimation,demand.estimation)
  
  # Counterfactual GSE costs?
  if(!is.null(RGSE)) { 
    supply.estimation[,MC_NONLIN_GSE := MC_NONLIN_GSE + RGSE]
  }
  
  ECAP_THRESHHOLD = 0.01
  
  # Get marginal costs for banks
  bank.data       <- fread('../data/final/capitalization_data_sep9.csv')
  merged.bank <- merge(bank.data,supply.estimation[LENDER == 'b',c('YEAR','LENDER','JUMBO','FICO','PURPOSE','MC_LIN_YEAR','MC_LIN_RISK_JUMBO','MC_LIN_LABOR','MC_NONLIN_GSE','MC_NONLIN_RE','MC_NONLIN_PSI'),with=F],by=c('YEAR','JUMBO','PURPOSE'),allow.cartesian = T)
  merged.bank[,ACTIVE_BASELINE := 1]
  
  if(!is.null(CR)) {
    merged.bank[,rho_hat := CR]
  }
  
  if(ALLOW_BANK_EXIT) {
    merged.bank[,ACTIVE_CF := (BANK_CR - rho_hat) > ECAP_THRESHHOLD]
  } else {
    merged.bank[,ACTIVE_CF := 1]
  }
  
  merged.bank[,ECAP_BANK := pmax(ECAP_THRESHHOLD,BANK_CR - rho_hat)]
  merged.bank[,mc_jumbo_nonlin_balance      := pmin(BANK_CR^2 * MC_NONLIN_RE * MC_NONLIN_PSI * xi_j * (ECAP_BANK)^(-1 * (MC_NONLIN_PSI + 1)),20)] # Just don't let anyone go crazy---they'd drop out anyway...
  merged.bank[,mc_conforming_nonlin_balance := BANK_CR^2 * MC_NONLIN_RE * MC_NONLIN_PSI * xi_g * (ECAP_BANK)^(-1 * (MC_NONLIN_PSI + 1))]
  merged.bank[,mc_conforming_nonlin_gse          := MC_NONLIN_GSE]
  merged.bank[JUMBO == 'j',mc_nonlin := pmin(mc_jumbo_nonlin_balance,mc_conforming_nonlin_gse + JUMBO_SEC_RATE,mc_conforming_nonlin_gse + EQUITY_ISSUANCE_RATE)] 
  merged.bank[JUMBO == 'c',mc_nonlin := pmin(mc_conforming_nonlin_balance,mc_conforming_nonlin_gse)]
  merged.bank[JUMBO == 'j',held := 1] # If you issue equity you're still holding it on balance sheet.....
  merged.bank[JUMBO == 'j',alternative := 0]
  merged.bank[JUMBO == 'j' & (mc_conforming_nonlin_gse + JUMBO_SEC_RATE) < mc_jumbo_nonlin_balance,alternative := 1]
  merged.bank[JUMBO == 'j' & (mc_conforming_nonlin_gse + EQUITY_ISSUANCE_RATE) < mc_jumbo_nonlin_balance,alternative := 1]
  merged.bank[JUMBO == 'j' & (mc_conforming_nonlin_gse + JUMBO_SEC_RATE) < mc_jumbo_nonlin_balance,held := 0] # Jumbo securitization
  merged.bank[JUMBO == 'c',held := 0]
  merged.bank[JUMBO == 'c' & mc_conforming_nonlin_balance < mc_conforming_nonlin_gse,held := 1]
  merged.bank[,MC := MC_LIN_YEAR + MC_LIN_RISK_JUMBO + MC_LIN_LABOR + mc_nonlin]
  bank.mc <- merged.bank[,j=list(LENDER = 'b',PCT_REMAIN = weighted.mean(ACTIVE_CF,w=N_ORIG) / weighted.mean(ACTIVE_BASELINE,w=N_ORIG),
                                 HELD = weighted.mean(held,w=N_ORIG * ACTIVE_CF),
                                 ALTERNATIVE = weighted.mean(alternative,w=N_ORIG * ACTIVE_CF),
                                 MC   = weighted.mean(MC,  w=N_ORIG * ACTIVE_CF)),by=c('YEAR','CBSA','PURPOSE','FICO','JUMBO')]
  
  # Combine it with fintech and non-fintech
  #non-bank
  nb.supply <- supply.estimation[LENDER != 'b']
  nb.supply[,HELD := 0]
  nb.supply[,PCT_REMAIN := 1]
  nb.supply[,ALTERNATIVE := 0]
  nb.supply[,MC := MC_LIN_LABOR + MC_LIN_YEAR + MC_NONLIN_GSE + MC_LIN_RISK_JUMBO]
  nb <- merge(estimation.data[LENDER != 'b',c('YEAR','CBSA','PURPOSE','FICO','JUMBO','LENDER')],nb.supply[,c('YEAR','FICO','PURPOSE','LENDER','HELD','ALTERNATIVE','PCT_REMAIN','MC')],by=c('YEAR','FICO','LENDER','PURPOSE'),allow.cartesian = T)
  
  combined <- rbind(nb,bank.mc)
  combined <- combined[order(YEAR,CBSA,PURPOSE,FICO,LENDER,JUMBO)]
  combined[is.na(MC),MC := 100]
  combined[is.na(PCT_REMAIN), PCT_REMAIN := 0]
  combined[PCT_REMAIN == 0,PCT_REMAIN := .01]
  
  capitalization.data <- combined
  
  supply.data <- merge(capitalization.data,market.data[,c('YEAR','CBSA','PURPOSE','FICO','MARKET_ID','JUMBO','LENDER','j','N_LENDERS')],by=c('YEAR','CBSA','PURPOSE','FICO','JUMBO','LENDER'))
  supply.data[,N_LENDERS := N_LENDERS * PCT_REMAIN]
  supply.data <- supply.data[,c('MARKET_ID','YEAR','CBSA','FICO','PURPOSE','j','LENDER','N_LENDERS','JUMBO','MC','HELD','ALTERNATIVE'),with=F]
  
  
  demand.data <- market.data[,c('MARKET_ID','YEAR','CBSA','PURPOSE','FICO','market.log_income.mean','market.log_income.std',
                                'market.log_price.mean','market.log_price.std','market.conforming_loan_limit','market.income_price.corr',
                                'market.log_income.mean.BAR','market.log_price.mean.BAR','MARKET_SIZE','MARKET_SHARE_DATA','RATE','market.log_loan_size.mean','LENDER','JUMBO','j','RATE_INST','delta','xi'),with=F]

  merged.cf.data <- merge(supply.data,demand.data,by=c('MARKET_ID','YEAR','CBSA','PURPOSE','FICO','j','LENDER','JUMBO'))
  
  
  # If Jumbo securitization...
  if(JUMBO_SEC_RATE < 999) {
    nb.rows <- merged.cf.data[LENDER != 'b']
    nb.rows[,j := j * 100]
    nb.rows[,JUMBO := 'j']
    nb.rows[,ALTERNATIVE := 1]
    nb.rows[,MC := MC + JUMBO_SEC_RATE] # Idea is you can securitize under exactly the same process but at some spread over what you'd do if GSE.
    merged.cf.data <- rbind(merged.cf.data,nb.rows)
    merged.cf.data <- merged.cf.data[order(MARKET_ID,j)]
    
  }
  
  toReturn <- list(market.data = merged.cf.data,linear.model = linear.model,demand.pars = demand.pars,raw.draws = raw.draws)
  
  return(toReturn)
}


estimation = 'estimation_final'

CR = NULL
LIM = NULL
RGSE = NULL
NAME = 'baseline'

runBootstraps(estimation)
